The script of the function spat_temp_index can be found in the following repository and supplementary material

Link to the functions spat_temp_index: Link to supplementary material:

# Load functions to run the all the simulations
setwd("C:/Users/David CM/Dropbox/DAVID DOC/LLAM al DIA/1. FEHM coses al DIA/4. Mecodispers Spa-Tem/MECODISPER SpaTem")
source("SpaTemp_HOBOS_function.R")

1. Charging & depurating HOBOS Dataset

The initial step is to charge the datasets containing the names of each monitored site containing the stream name, the ID code and the latitude and longitude. This dataset will be used to localize each one of the data loggers along the river (the “Longlat_Rius.csv” can be found in the repository - LINK)

The following step is to upload the HOBOS_sites information, which corresponds to the A) Water permanence database (Figure 1). We upload individually the tables for each stream (one sheet for each stream). These “.csv” documents can be found in the online repository - LINK and contain the information about each daily wet/dry status of each stream.

IMPORTANT INFORMATION: For each one of these “.csv” documents the data loggers must be ordered in Upstream to Downstream direction and must have the following structure:

data.frame(
MonitoredDays= c("Day1","Day2","Day3","Day4"), 
UpstreamSite1=c(0,1,1,1), 
StreamSite2=c(0,1,1,1),StreamSite3=c(1,1,1,1),StreamSite4=c(1,1,1,1),
DownstreamSite5=c(1,1,1,1) 
)

Finally, we upload the “Stream_order” database. This database contains the Upstream - Downstream “order” of each datalogger (1 being Upstream). The file “Stream Order.csv” can be found in the online repository - LINK. Later, we check and correct some coordinates to make sure that each datalogger has a unique coordinate and plot them individually to assess if they are correctly placed in space.

With this initial process we will depurate and charge the two most important informations required for the function: 1. HOBOS_dataset 2. Sites_coordinates

Sites <- read.csv("Raw_HOBOS_Database/Longlat_Rius.csv", header = T, sep = ";", dec = ",")
colnames(Sites) <- c("Riera", "Codi_HOBO","Latitud","Longitud")

#Correction for matching HOBOS -- collecting the  coordinates of the HOBOS that we have  
length(which(as.matrix(dist(Sites[,3:4]))==0))
length(diag(as.matrix(dist(Sites[,3:4]))))
plot(Sites$Longitud, Sites$Latitud)

# We upload the HOBO dataset (1 and 0 defining wet/dry moments)
# We upload the order of the rivers from Upstream to Downstream to arrange the HOBOS into the desired order
## IMPORTANT: This order will also define the direction of the river! 
## The first HOBO in the row of the "Sites_list" must be the first HOBO in the Column of "HOBOS_sites"
HOBOS_sites <- list(
  read.csv("Raw_HOBOS_Database/CA_HOBOS_data.csv", header = T, sep = ";"),
  read.csv("Raw_HOBOS_Database/M_HOBOS_data.csv", header = T, sep = ";"),
  read.csv("Raw_HOBOS_Database/R_HOBOS_data.csv", header = T, sep = ";"),
  read.csv("Raw_HOBOS_Database/SA_HOBOS_data.csv", header = T, sep = ";"),
  read.csv("Raw_HOBOS_Database/SC_HOBOS_data.csv", header = T, sep = ";"),
  read.csv("Raw_HOBOS_Database/T_HOBOS_data.csv", header = T, sep = ";"),
  read.csv("Raw_HOBOS_Database/VH_HOBOS_data.csv", header = T, sep = ";"))

#Stream order to check
Stream_order <- read.csv("Raw_HOBOS_Database/Stream Order.csv", header = T, sep = ";")

library(dplyr)
# PLotting HOBOS altogether (make the plot window bigger)
par(mfrow=c(3,3))
local_hobos_list <- list()
Sites_list <- list()
for (sit in 1:length(HOBOS_sites)) {
  colnames(HOBOS_sites[[sit]])[1] <- c("Day")
  names_hobos <- colnames(HOBOS_sites[[sit]])[2:length(colnames(HOBOS_sites[[sit]]))]
  local_hobos <- c()
  Sites_list[[sit]] <- Sites%>%filter(Codi_HOBO%in%names_hobos)%>%
                               left_join(Stream_order, by="Codi_HOBO")%>%
                               arrange(UtoD)%>%
                               dplyr::select(Riera,Codi_HOBO,Latitud,Longitud)
plot(Sites_list[[sit]]$Longitud, Sites_list[[sit]]$Latitud)
}
par(mfrow=c(1,1))

2. Dates selection

To select different “time windows” in order to calculate ST indicies for the the whole monitored period (513 days) or until the Autumn sampling (18-21/11/2019; 480 days) we upload the biological information files and retrieve and modify the dates in order to easily select them and calculate the difference between the beginning (“date_HOBOS” = 2018-07-26) and the other two dates, the end (“date_HOBOS_end”=2019-12-20) and the autumn sampling (“date_correct_Autumn”).

As each row of the matrix “HOBOS_sites” corresponds to a unique day, the difference (in days) between two dates (beginning vs end or beginning vs Autumn sampling) will deliver the number of rows that we have to select from the “HOBOS_sites” in order to account for the desired “time-window”. Therefore, after this second step we will have the correct HOBOS_sites matrix corresponding to the time-window that we want to focus.

# We upload the biological data just to see the Date of the interval that we want to capture 
### Autumn dataset - November 2019
BDD <- read.csv2("BiolData/Matriz_otoño.csv", sep=";")
date_correct_Autumn <- min(unique(BDD$Date))
# Edit the dates according the desired format to be entered in the following steps
date_correct_Autumn <- strsplit(date_correct_Autumn,"/") # Split the numbers by the "/"
date_correct_Autumn <- paste(date_correct_Autumn[[1]][3],"-", # Paste the numbers in the correct order and separated by "-" 
                             date_correct_Autumn[[1]][2],"-",
                             date_correct_Autumn[[1]][1],sep = "")

date_correct_Autumn
## All hobos begun at the same date
date_HOBOS <- "2018-07-26"

## All hobos finish at the same date
date_HOBOS_end <- "2019-12-20"

#Beginning
bd <- as.Date(date_HOBOS)
#End
ed <- as.Date(date_HOBOS_end)

# Difference in number of days will correspond to the number of rows to be selected. 
time_window_beg <- as.numeric(difftime(bd+1, date_HOBOS, units = "days")) 
time_window_beg <- round(time_window_beg)
time_window_end <- as.numeric(difftime(ed, date_HOBOS, units = "days"))
time_window_end <- round(time_window_end)

# We select the corresponding number of days for each stream
for (site in 1:length(HOBOS_sites)) {
  HOBOS_sites[[site]] <- HOBOS_sites[[site]][time_window_beg:time_window_end,]
}

3. Distance matric creation

Finally, if we wish to calculate weighted indices we have to create a weighted pairwise matrix with the desired distance. In our case we caluclate the euclidean distances based on the longitude and latitude. We use the coordinates information contained in “Sites_list” object.

By doing this we obtain the last required information for running the function, the dist_matrices, here named “eucl_dist_matrices” as we used the direct euclidean distances between monitored sites.

# Must be a list of distances matrix for each river. 
# Must be a symmetric matrix (upper and lower triangles must be equal). The function already selects which distances 
#are selected in each case.  
eucl_dist_matrices <- list()
for (river in 1:length(Sites_list)) {
  eucl_dist_matrices[[river]] <-as.matrix(dist(Sites_list[[river]][,4:3],diag = T,upper = T)) 
}

Scenario A1 - Binnary directed stream

# Building a directed non weighted network
# Value of link=1
# Value of NO link=0
Dir_NonW_Net <- spat_temp_index(HOBOS_dataset=HOBOS_sites,
                                Sites_coordinates=Sites_list,
                                direction="directed", 
                                weighting=FALSE,
                                value_S_LINK=1,
                                value_T_LINK=1,
                                value_NO_S_link=0,
                                value_NO_T_link=0,
                                Network_variables=TRUE,
                                print.plots=TRUE,
                                print.directory="Figure/")

Scenario A2 - Weighted directed stream

# Building a directed Weighted network
# Value of link= 0.1
# Value of NO link=1
# - I am evaluationg "resistance" to move from one to another. 
# SO, higher numbers mean High resistance
Dir_WEIG_Net <- spat_temp_index(HOBOS_dataset=HOBOS_sites,
                                Sites_coordinates=Sites_list,
                                direction="directed", 
                                weighting=TRUE,
                                dist_matrices = eucl_dist_matrices,
                                value_S_LINK=0.1,
                                value_T_LINK=0.1,
                                value_NO_S_link=1,
                                value_NO_T_link=1,
                                Network_variables=TRUE,
                                print.plots=TRUE,
                                print.directory="Figure/")

Scenario B1 - Binnary undirected stream

# Building a Undirected non weighted network
# Value of link=1
# Value of NO link=0
UnD_NonW_Net <- spat_temp_index(HOBOS_dataset=HOBOS_sites,
                                Sites_coordinates=Sites_list,
                                direction="undirected", 
                                weighting=FALSE,
                                value_S_LINK=1,
                                value_T_LINK=1,
                                value_NO_S_link=0,
                                value_NO_T_link=0,
                                Network_variables=TRUE,
                                print.plots=TRUE,
                                print.directory="Figure/")

Scenario B2 - Weighted undirected stream

# Building a Undirected Weighted network
# Value of link= 0.1
# Value of NO link=1
# - I am evaluationg "resistance" to move from one to another. 
# SO, higher numbers mean High resistance
UnD_WEIG_Net <- spat_temp_index(HOBOS_dataset=HOBOS_sites,
                                Sites_coordinates=Sites_list,
                                direction="undirected", 
                                weighting=TRUE,
                                dist_matrices = eucl_dist_matrices,
                                value_S_LINK=0.1,
                                value_T_LINK=0.1,
                                value_NO_S_link=1,
                                value_NO_T_link=1,
                                Network_variables=TRUE,
                                print.plots=TRUE,
                                print.directory="Figure/")

A1 Scenario - Binnary directed stream

# Directed NONWEIGHTED river network OUTPUTS ####
# Spatiotemporal matrix 
NonW_ST_matrix_out_out <- Dir_NonW_Net$Dir_NonW_ST_matrix

# Spatiotemporal connectivity 
NonW_ST_connectivity_value <- Dir_NonW_Net$Dir_NonW_ST_con

# Spatiotemporal Out.closenness 
NonW_ST_directed_Ocloseness_rivers <- Dir_NonW_Net$Dir_NonW_ST_Ocl
NonW_ST_directed_Oclo_mean <- list()
NonW_ST_directed_Oclo_sd <- list()
for (river in 1:length(NonW_ST_directed_Ocloseness_rivers)) {
  NonW_ST_directed_Oclo_mean[[river]] <- apply(NonW_ST_directed_Ocloseness_rivers[[river]],2,mean)
  NonW_ST_directed_Oclo_sd[[river]] <- apply(NonW_ST_directed_Ocloseness_rivers[[river]],2,sd)
}

#NonW_ST_oclo_plot
# Spatiotemporal ALL.closenness 
NonW_ST_directed_Allcloseness_rivers <- Dir_NonW_Net$Dir_NonW_ST_Acl
NonW_ST_directed_Allclo_mean <- list()
NonW_ST_directed_Allclo_sd <- list()
for (river in 1:length(NonW_ST_directed_Allcloseness_rivers)) {
  NonW_ST_directed_Allclo_mean[[river]] <- apply(NonW_ST_directed_Allcloseness_rivers[[river]],2,mean)
  NonW_ST_directed_Allclo_sd[[river]] <- apply(NonW_ST_directed_Allcloseness_rivers[[river]],2,sd)
}

#NonW_ST_Allclo_plot
# Spatiotemporal Betweenness 
NonW_ST_directed_betweennes_rivers <- Dir_NonW_Net$Dir_NonW_ST_Bet
NonW_ST_directed_betw_mean <- list()
NonW_ST_directed_betw_sd <- list()
for (river in 1:length(NonW_ST_directed_betweennes_rivers)) {
  NonW_ST_directed_betw_mean[[river]] <- apply(NonW_ST_directed_betweennes_rivers[[river]],2,mean)
  NonW_ST_directed_betw_sd[[river]] <- apply(NonW_ST_directed_betweennes_rivers[[river]],2,sd)
}

#NonW_ST_betw_plot
NonW_ST_directed_output <- data.frame(NonW_Dir_con=unlist(NonW_ST_connectivity_value),
                                      NonW_Dir_Ocl=unlist(NonW_ST_directed_Oclo_mean),
                                      NonW_Dir_Acl=unlist(NonW_ST_directed_Allclo_mean),
                                      NonW_Dir_Bet=unlist(NonW_ST_directed_betw_mean))

A2 Scenario - Weighted directed stream

# Directed WEIGHTED river network OUTPUTS ####
# Spatiotemporal matrix 
WEIG_ST_matrix_out_out <- Dir_WEIG_Net$Dir_WEIG_ST_matrix

# Spatiotemporal connectivity 
WEIG_ST_connectivity_value <- Dir_WEIG_Net$Dir_WEIG_ST_con

# Spatiotemporal Out.closenness 
WEIG_ST_directed_Ocloseness_rivers <- Dir_WEIG_Net$Dir_WEIG_ST_Ocl
WEIG_ST_Oclo_mean <- list()
WEIG_ST_Oclo_sd <- list()
for (river in 1:length(WEIG_ST_directed_Ocloseness_rivers)) {
  WEIG_ST_Oclo_mean[[river]] <- apply(WEIG_ST_directed_Ocloseness_rivers[[river]],2,mean)
  WEIG_ST_Oclo_sd[[river]] <- apply(WEIG_ST_directed_Ocloseness_rivers[[river]],2,sd)
}

# Spatiotemporal All.closenness 
WEIG_ST_directed_Allcloseness_rivers <- Dir_WEIG_Net$Dir_WEIG_ST_Acl
WEIG_ST_Allclo_mean <- list()
WEIG_ST_Allclo_sd <- list()
for (river in 1:length(WEIG_ST_directed_Allcloseness_rivers)) {
  WEIG_ST_Allclo_mean[[river]] <- apply(WEIG_ST_directed_Allcloseness_rivers[[river]],2,mean)
  WEIG_ST_Allclo_sd[[river]] <- apply(WEIG_ST_directed_Allcloseness_rivers[[river]],2,sd)
}

# Spatiotemporal Betweenness 
WEIG_ST_directed_betweennes_rivers <- Dir_WEIG_Net$Dir_WEIG_ST_Bet
WEIG_ST_betw_mean <- list()
WEIG_ST_betw_sd <- list()
for (river in 1:length(WEIG_ST_directed_betweennes_rivers)) {
  WEIG_ST_betw_mean[[river]] <- apply(WEIG_ST_directed_betweennes_rivers[[river]],2,mean)
  WEIG_ST_betw_sd[[river]] <- apply(WEIG_ST_directed_betweennes_rivers[[river]],2,sd)
}

WEIG_ST_directed_output <- data.frame(WEIG_Dir_con=unlist(WEIG_ST_connectivity_value),
                                      WEIG_Dir_Ocl=unlist(WEIG_ST_Oclo_mean),
                                      WEIG_Dir_Acl=unlist(WEIG_ST_Allclo_mean),
                                      WEIG_Dir_Bet=unlist(WEIG_ST_betw_mean))

Scenario B1 - Binnary undirected stream

# Undirected river network OUTPUTS ####
# Spatiotemporal matrix 
Un_NonW_ST_matrix_out_out <- UnD_NonW_Net$UnD_NonW_ST_matrix

# Spatiotemporal matrix 
Un_NonW_ST_connectivity_value <- UnD_NonW_Net$UnD_NonW_ST_con

# Spatiotemporal All closenness 
Un_NonW_ST_directed_Ocloseness_rivers <- UnD_NonW_Net$UnD_NonW_ST_Ocl
Un_NonW_ST_Oclo_mean <- list()
Un_NonW_ST_Oclo_sd <- list()
for (river in 1:length(Un_NonW_ST_directed_Ocloseness_rivers)) {
  Un_NonW_ST_Oclo_mean[[river]] <- apply(Un_NonW_ST_directed_Ocloseness_rivers[[river]],2,mean)
  Un_NonW_ST_Oclo_sd[[river]] <- apply(Un_NonW_ST_directed_Ocloseness_rivers[[river]],2,sd)
}

# Spatiotemporal Out closenness 
Un_NonW_ST_directed_Allcloseness_rivers <- UnD_NonW_Net$UnD_NonW_ST_Acl
Un_NonW_ST_Allclo_mean <- list()
Un_NonW_ST_Allclo_sd <- list()
for (river in 1:length(Un_NonW_ST_directed_Allcloseness_rivers)) {
  Un_NonW_ST_Allclo_mean[[river]] <- apply(Un_NonW_ST_directed_Allcloseness_rivers[[river]],2,mean)
  Un_NonW_ST_Allclo_sd[[river]] <- apply(Un_NonW_ST_directed_Allcloseness_rivers[[river]],2,sd)
}

# Spatiotemporal Betweenness 
Un_NonW_ST_directed_betweennes_rivers <- UnD_NonW_Net$UnD_NonW_ST_Bet
Un_NonW_ST_betw_mean <- list()
Un_NonW_ST_betw_sd <- list()
for (river in 1:length(Un_NonW_ST_directed_betweennes_rivers)) {
  Un_NonW_ST_betw_mean[[river]] <- apply(Un_NonW_ST_directed_betweennes_rivers[[river]],2,mean)
  Un_NonW_ST_betw_sd[[river]] <- apply(Un_NonW_ST_directed_betweennes_rivers[[river]],2,sd)
}

Un_NonW_ST_output <- data.frame(Un_NonW_con=unlist(Un_NonW_ST_connectivity_value),
                                Un_NonW_Ocl=unlist(Un_NonW_ST_Oclo_mean),
                                Un_NonW_Acl=unlist(Un_NonW_ST_Allclo_mean),
                                Un_NonW_Bet=unlist(Un_NonW_ST_betw_mean))

Scenario B2 - Weighted undirected stream

# Undirected WEIGHTED river network OUTPUTS ####
# Spatiotemporal matrix 
Un_WEIG_ST_matrix_out_out <- UnD_WEIG_Net$UnD_WEIG_ST_matrix

# Spatiotemporal matrix 
Un_WEIG_ST_connectivity_value <- UnD_WEIG_Net$UnD_WEIG_ST_con

# Spatiotemporal All closenness 
Un_WEIG_ST_directed_Ocloseness_rivers <- UnD_WEIG_Net$UnD_WEIG_ST_Ocl
Un_WEIG_ST_Oclo_mean <- list()
Un_WEIG_ST_Oclo_sd <- list()
for (river in 1:length(Un_WEIG_ST_directed_Ocloseness_rivers)) {
  Un_WEIG_ST_Oclo_mean[[river]] <- apply(Un_WEIG_ST_directed_Ocloseness_rivers[[river]],2,mean)
  Un_WEIG_ST_Oclo_sd[[river]] <- apply(Un_WEIG_ST_directed_Ocloseness_rivers[[river]],2,sd)
}


# Spatiotemporal Out closenness 
Un_WEIG_ST_directed_Allcloseness_rivers <- UnD_WEIG_Net$UnD_WEIG_ST_Acl
Un_WEIG_ST_Allclo_mean <- list()
Un_WEIG_ST_Allclo_sd <- list()
for (river in 1:length(Un_WEIG_ST_directed_Allcloseness_rivers)) {
  Un_WEIG_ST_Allclo_mean[[river]] <- apply(Un_WEIG_ST_directed_Allcloseness_rivers[[river]],2,mean)
  Un_WEIG_ST_Allclo_sd[[river]] <- apply(Un_WEIG_ST_directed_Allcloseness_rivers[[river]],2,sd)
}

# Spatiotemporal Betweenness 
Un_WEIG_ST_directed_betweennes_rivers <- UnD_WEIG_Net$UnD_WEIG_ST_Bet
Un_WEIG_ST_Betclo_mean <- list()
Un_WEIG_ST_Betclo_sd <- list()
for (river in 1:length(Un_WEIG_ST_directed_betweennes_rivers)) {
  Un_WEIG_ST_Betclo_mean[[river]] <- apply(Un_WEIG_ST_directed_betweennes_rivers[[river]],2,mean)
  Un_WEIG_ST_Betclo_sd[[river]] <- apply(Un_WEIG_ST_directed_betweennes_rivers[[river]],2,sd)
}

Un_WEIG_ST_output <- data.frame(Un_WEIG_con=unlist(Un_WEIG_ST_connectivity_value),
                                Un_WEIG_Ocl=unlist(Un_WEIG_ST_Oclo_mean),
                                Un_WEIG_Acl=unlist(Un_WEIG_ST_Allclo_mean),
                                Un_WEIG_Bet=unlist(Un_WEIG_ST_Betclo_mean))

# We add the stream ID and the order (UtoD) to each output to later analyse them. 
# We retrieve this information from the "Stream_order" dataset where each HOBO code, Stram and its position on the 
#upstream downstream direction were registered
# This two following vectors are created for being used in other scripts
HOB_riv_ID<- Stream_order$ï..Riera
ups_dos <- Stream_order$UtoD

NonW_ST_directed_out <- data.frame(ID=Stream_order$ï..Riera,DtoU=Stream_order$UtoD,
                                   NonW_ST_directed_output)%>%
                        mutate(ID_UpDo=paste(ID,"_",DtoU,sep = ""))

WEIG_ST_directed_out <- data.frame(ID=Stream_order$ï..Riera,DtoU=Stream_order$UtoD,
                                   WEIG_ST_directed_output)%>%
                        mutate(ID_UpDo=paste(ID,"_",DtoU,sep = ""))

Un_NonW_ST_out <- data.frame(ID=Stream_order$ï..Riera,DtoU=Stream_order$UtoD,
                             Un_NonW_ST_output)%>%
                        mutate(ID_UpDo=paste(ID,"_",DtoU,sep = ""))

Un_WEIG_ST_out <- data.frame(ID=Stream_order$ï..Riera,DtoU=Stream_order$UtoD,
                             Un_WEIG_ST_output)%>%
                        mutate(ID_UpDo=paste(ID,"_",DtoU,sep = ""))


# Final values HOBOS dataframe
NonW_ST_matrix_out_out
WEIG_ST_matrix_out_out
Un_NonW_ST_matrix_out_out
Un_WEIG_ST_matrix_out_out

NonW_ST_directed_out
WEIG_ST_directed_out
Un_NonW_ST_out
Un_WEIG_ST_out

Initial and required data to proceed (previously crated)

# Cunillera_palette - Available at https://github.com/Cunillera-Montcusi/Cunillera_Pallete
source("C:/Users/David CM/Dropbox/DAVID DOC/LLAM al DIA/CUNILLERA_palette.R")
setwd("C:/Users/David CM/Dropbox/DAVID DOC/LLAM al DIA/1. FEHM coses al DIA/4. Mecodispers Spa-Tem/MECODISPER SpaTem")

## IMPORTANT MESSAGE: YOU NEED TO RUN the treatment script
# Needed values from the script "SpaTemp_HOBOS_treatment.R"
## RUN BEFORE PROCEEDING WITH THIS ONE

## You need to run the previous script in order to built the HOBOS dataset (a list with the HOBOS info for each river)
HOBOS_sites
# This two other objects are also built in the "SpaTemp_HOBOS_treatment.R" script. They contain the names 
## of each river and the order upstream to downstream within it. 
HOB_riv_ID
ups_dos

# Matrix ST outputs obtained from "SpaTemp_HOBOS_treatment.R"
NonW_ST_matrix_out_out
WEIG_ST_matrix_out_out
Un_NonW_ST_matrix_out_out
Un_WEIG_ST_matrix_out_out

# INDIVIDUAL ST outputs obtained from "SpaTemp_HOBOS_treatment.R"
NonW_ST_directed_out
WEIG_ST_directed_out
Un_NonW_ST_out
Un_WEIG_ST_out

Individual HOBO values and comparisons with old metrics - STcon tratment

We calculate a different set of metrics that where previously used for this type of dataset to evaluate drought patterns. For each site we calculated the total duration of drying events (i.e. the total number of dry days; hereafter TotDur), the frequency of drying events (i.e. the number of changes from wet to dry days; hereafter TotNum), and the average length of each drying event (the average number of dry days for each drying event; hereafter TotLeng).

The script to calculate these metrics can be found in the “Old_HOBOS_calculation.R” file found in the following repository - LINK

source("C:/Users/David CM/Dropbox/DAVID DOC/LLAM al DIA/1. FEHM coses al DIA/4. Mecodispers Spa-Tem/MECODISPER SpaTem/Old_HOBOS_calculation.R")
Old_HOBOS_comb

Once obtained, we compare them with the STcon indices obtained for the four scenarios studied (PCA analysis-Figure 3) and analyse the correlation among them (Figure S3 ). More information in the article DOI:

Figure 3: PCA plot for STcon indices for each one of the four scenarios: Directed binary (STcon.A1), Directed Weighted (STcon.A2), Undirected Binary (STcon.B1), Undirected Weighted (STcon.B2) and the previously existing indices for water permanence data: the total duration of drying events (TotDur), the frequency of drying events (TotNum), and the average length of each drying event (TotLeng). All the data is calculated based on the total monitored time window monitored (513 days).

# Figure 3 
library(vegan)
library(ggfortify)
data_PCA <- data.frame(data.frame("STcon A1"=NonW_ST_directed_out$NonW_Dir_con,
                                  "STcon A2"=WEIG_ST_directed_out$WEIG_Dir_con,
                                  "STcon B1"=Un_NonW_ST_out$Un_NonW_con, 
                                  "STcon B2"=Un_WEIG_ST_out$Un_WEIG_con),
            Old_HOBOS_comb[,3:5])

PCA_info <- prcomp(data_PCA,center = T, scale. = T)
a <- data.frame(
     adonis2(PCA_info$x[,2]~NonW_ST_directed_out$ID*NonW_ST_directed_out$DtoU,
        method = "euclidean")
      )
rownames(a) <- c("Stream ID","Upstr-Downstr", "Interaction","Residual","Total")
caption_text <- paste(rownames(a)[1],"pvalue=",a[1,5]," ",
                       rownames(a)[2],"pvalue=",a[2,5]," ",
                       rownames(a)[3],"pvalue=",a[3,5], sep = " ")
autoplot(prcomp(data_PCA,
  center = T, scale. = T),
  loadings=T, loadings.label = TRUE, loadings.label.size = 6, shape = F, label=F,
  loadings.colour = 'grey15', loadings.label.colour="black")+
  geom_point(alpha=0.5,shape=21,
             color="grey30", 
             aes(size=NonW_ST_directed_out$DtoU ,fill=NonW_ST_directed_out$ID))+
  stat_ellipse(aes(colour=NonW_ST_directed_out$ID),type = "t", linetype=2,size=1, alpha=0.5)+
  scale_fill_CUNILLERA(palette = "LGTBI", name="Stream ID")+
  scale_color_CUNILLERA(palette = "LGTBI", name="Stream ID")+
  scale_size(name="Upstr-Downstr")+
  scale_x_continuous(limits = c(-0.35,0.4))+
  scale_y_continuous(limits = c(-0.32,0.32))+
  labs(caption = caption_text)+
  theme_bw()

Figure S3: Correlation plot among all the indices obtained for the four scenarios and the previously used drought-related metrics.

# Figure 3S
corrmorant::corrmorant(cbind("STcon 1.1."=NonW_ST_directed_out$NonW_Dir_con,
                             "log(STcon) 1.2."=log(WEIG_ST_directed_out$WEIG_Dir_con+1),
                             "STcon 2.1"=Un_NonW_ST_out$Un_NonW_con,
                             "log(STcon) 2.2"=log(Un_WEIG_ST_out$Un_WEIG_con+1),
                             "log(TotDur)"= log(Old_HOBOS_comb[,3]+1),
                             "log(TotNum)"= log(Old_HOBOS_comb[,4]+1),
                             "log(TotLeng)"= log(Old_HOBOS_comb[,5]+1)),
                            rescale = "by_sd",corr_method = "spearman")

Stream characterization using STcon

The obtained indices can be used to characterise the river spatiotemporally by ploting their STcon vlaues along the stream direction (Figure S4) for example, wich allows to detect more disconnected sites or the impact of determined sites on the rest of the network (e.g. directed networks).

For these characterization one can use either the STcon values, Out closennes, All closennes or Betweenness depending on the main aim or the specific feature that has to be targeted.

library(gridExtra)
four_approxim <- list(NonW_ST_directed_out, 
                      WEIG_ST_directed_out,
                      Un_NonW_ST_out, 
                      Un_WEIG_ST_out)

names_approxim <- c("Dir_NonW_ST", 
                    "Dir_WEIG_ST",
                    "Un_NonW_ST", 
                    "Un_WEIG_ST")

for (a in 1:length(four_approxim)) {
rows_filter <- unlist(c(four_approxim[[a]]%>%
                          mutate(files=1:nrow(.))%>%
                          group_by(ID)%>%
                          mutate(maxim=max(DtoU), value=DtoU/maxim)%>%
                          filter(value==1)%>%
                          ungroup()%>%
                          dplyr::select(files)))
    
grid.arrange(
  ggplot(data = four_approxim[[a]][-rows_filter,])+
    geom_point(aes(x=DtoU,y=four_approxim[[a]][-rows_filter,3],fill=ID),color="grey30", alpha=0.5, shape=21, size=2)+
    geom_line(aes(x=DtoU,y=four_approxim[[a]][-rows_filter,3],color=ID),alpha=0.3, linetype=2)+
    geom_smooth(aes(x=DtoU,y=four_approxim[[a]][-rows_filter,3],color=ID),alpha=0.2, method = "lm",se = F)+
    scale_fill_CUNILLERA(palette = "LGTBI")+scale_color_CUNILLERA(palette = "LGTBI")+
    xlab("HOBO relative position Upstream-Downstream")+ylab("ST Connecticity")+theme_classic()+theme(legend.position="none"),
  
  ggplot(data = four_approxim[[a]])+
    geom_point(aes(x=DtoU,y=four_approxim[[a]][,4] ,fill=ID),color="grey30", alpha=0.5, shape=21, size=2)+
    geom_line(aes(x=DtoU,y=four_approxim[[a]][,4] ,color=ID),alpha=0.3, linetype=2)+
    geom_smooth(aes(x=DtoU,y=four_approxim[[a]][,4] ,color=ID),alpha=0.2, method = "lm",se = F)+
    scale_fill_CUNILLERA(palette = "LGTBI")+scale_color_CUNILLERA(palette = "LGTBI")+
    xlab("HOBO relative position Upstream-Downstream")+ylab("ST Out closeness")+theme_classic()+theme(legend.position="none"),
  
  ggplot(data = four_approxim[[a]])+
    geom_point(aes(x=DtoU,y=four_approxim[[a]][,5]  ,fill=ID),color="grey30", alpha=0.5, shape=21, size=2)+
    geom_line(aes(x=DtoU,y=four_approxim[[a]][,5]  ,color=ID),alpha=0.3, linetype=2)+
    geom_smooth(aes(x=DtoU,y=four_approxim[[a]][,5]  ,color=ID),alpha=0.2, method = "lm",se = F)+
    scale_fill_CUNILLERA(palette = "LGTBI")+scale_color_CUNILLERA(palette = "LGTBI")+xlab("HOBO relative position Upstream-Downstream")+
    ylab("ST All closeness")+theme_classic()+theme(legend.position="none"),
  
  ggplot(data = four_approxim[[a]])+
    geom_point(aes(x=DtoU,y=four_approxim[[a]][,6]  ,fill=ID),color="grey30", alpha=0.5, shape=21, size=2)+
    geom_line(aes(x=DtoU,y=four_approxim[[a]][,6]  ,color=ID),alpha=0.3, linetype=2)+
    geom_smooth(aes(x=DtoU,y=four_approxim[[a]][,6]  ,color=ID),alpha=0.2, method = "lm",se = F)+
    scale_fill_CUNILLERA(palette = "LGTBI")+scale_color_CUNILLERA(palette = "LGTBI")+
    xlab("HOBO relative position Upstream-Downstream")+ylab("ST Betwenness")+theme_classic()+theme(legend.position="none"),
  
  get_legend(ggplot(NonW_ST_directed_out)+
               geom_point(aes(x=DtoU,y=four_approxim[[a]][,3],fill=ID),shape=21, size=2)+
               scale_fill_CUNILLERA(palette = "LGTBI", name="Stream ID")+theme_classic()+theme(legend.direction = "horizontal",legend.box="vertical")),
  nrow=3 ,ncol=2 ,top=names_approxim[a])
}

Matrix HOBOS values NMDS building - STconmat and stream characterization using both indices

For each scenario (A1, A2, B1 and B2) we calculate an NMDS using euclidean distances and ploting their values in order to characterize each one of the analysed streams in a bidimensional space.

Finally, we merge the obtained results in a unique figure where the two developed indices are shown. Figur 4 of the article (DOI: LINK) provides a complete image of streams spatiotemporal patterns based on both indices.

# Figure 4
  rows_filter <- unlist(c(four_approxim[[1]]%>%
                            mutate(files=1:nrow(.))%>%
                            group_by(ID)%>%
                            mutate(maxim=max(DtoU), value=DtoU/maxim)%>%
                            filter(value==1)%>%
                            ungroup()%>%
                            dplyr::select(files)))
  
  alpha_stream <- as.numeric(unlist(c(four_approxim[[1]]%>%
                  mutate(alpha_stream=ifelse(ID=="VH", 1, 0.3))%>%
                  ungroup()%>%
                  dplyr::select(alpha_stream))))
  alpha_HOBOS3 <- c(rep(0,58),1,rep(0,7))
  
  grid.arrange(
  #A1 line plot  
  arrangeGrob(
    ggplot(data = four_approxim[[1]][-rows_filter,])+
      geom_point(aes(x=DtoU,y=four_approxim[[1]][-rows_filter,3]),color="red",alpha=alpha_HOBOS3[-rows_filter], shape=16, size=5,)+
      geom_point(aes(x=DtoU,y=four_approxim[[1]][-rows_filter,3],fill=ID),color="grey30", alpha=alpha_stream[-rows_filter], shape=21, size=2)+
      geom_line(aes(x=DtoU,y=four_approxim[[1]][-rows_filter,3],color=ID),alpha=alpha_stream[-rows_filter], linetype=2)+
      geom_line(aes(x=DtoU,y=four_approxim[[1]][-rows_filter,3],color=ID, alpha=ID), 
                size=1.3,stat="smooth",method = "lm")+
      scale_alpha_manual(values = c(rep(0.3,6),0.9))+
      scale_fill_CUNILLERA(palette = "LGTBI")+
      scale_color_CUNILLERA(palette = "LGTBI")+
      labs(subtitle = "A1 - STcon", title = "A)")+
      xlab("HOBO relative position Upstream-Downstream")+
      ylab("ST Connecticity")+
      theme_classic()+theme(legend.position="none"),
    top = ""),
  
  #A1 NMDS plot  
  arrangeGrob(    
    ggplot(NonW_ST_directed_MatrixOut)+
      #geom_point(aes(x=x,y=y, fill=ID,size=DtoU),color="red",alpha=alpha_HOBOS3, shape=16, size=5,)+
      geom_point(aes(x=x, y=y, fill=ID,size=DtoU), shape=21,colour="grey10", alpha=alpha_stream)+
      stat_ellipse(aes(x=x, y=y, colour=ID),type = "t", linetype=2,size=1, alpha=0.5)+
      scale_fill_CUNILLERA(palette = "LGTBI")+
      scale_color_CUNILLERA(palette = "LGTBI")+
      labs(subtitle = "A1 - STconmat", title = "B)")+
      scale_y_continuous(label=NULL)+
      scale_x_continuous(label=NULL)+
      theme_classic()+theme(legend.position="none",axis.ticks = element_blank())+
      facet_grid(. ~ ID),
    top = ""),
    
  
  #A2 line plot  
  arrangeGrob( 
    ggplot(data = four_approxim[[2]][,])+
      geom_point(aes(x=DtoU,y=four_approxim[[2]][,3]),color="red",alpha=alpha_HOBOS3, shape=16, size=5,)+
      geom_point(aes(x=DtoU,y=four_approxim[[2]][,3],fill=ID),color="grey30", alpha=alpha_stream, shape=21, size=2)+
      geom_line(aes(x=DtoU,y=four_approxim[[2]][,3],color=ID),alpha=alpha_stream, linetype=2)+
      geom_line(aes(x=DtoU,y=four_approxim[[2]][,3],color=ID, alpha=ID), 
                size=1.3,stat="smooth",method = "lm")+
      scale_alpha_manual(values = c(rep(0.3,6),0.9))+
      scale_fill_CUNILLERA(palette = "LGTBI")+
      scale_color_CUNILLERA(palette = "LGTBI")+
      labs(subtitle = "A2 - STcon")+ 
      xlab("HOBO relative position Upstream-Downstream")+
      ylab("ST Connecticity")+
      scale_y_continuous(limits = c(0,15000))+
      theme_classic()+theme(legend.position="none"),
    top=""),
    
  #A2 NMDS plot
  arrangeGrob(
    ggplot(WEIG_ST_MatrixOut)+
      geom_point(aes(x=x, y=y, fill=ID,size=DtoU), shape=21,colour="grey10", alpha=alpha_stream)+
      stat_ellipse(aes(x=x, y=y, colour=ID),type = "t", linetype=2,size=1, alpha=0.5)+
      scale_fill_CUNILLERA(palette = "LGTBI")+
      scale_color_CUNILLERA(palette = "LGTBI")+
      scale_y_continuous(label=NULL)+
      scale_x_continuous(label=NULL)+
      labs(subtitle = "A2 - STconmat")+
      theme_classic()+theme(legend.position="none",axis.ticks = element_blank())+
    facet_grid(. ~ ID), top = ""),
  
  #B1 line plot
  arrangeGrob(
    ggplot(data = four_approxim[[3]][,])+
      geom_point(aes(x=DtoU,y=four_approxim[[3]][,3]),color="red",alpha=alpha_HOBOS3, shape=16, size=5,)+
      geom_point(aes(x=DtoU,y=four_approxim[[3]][,3],fill=ID),color="grey30", alpha=alpha_stream, shape=21, size=2)+
      geom_line(aes(x=DtoU,y=four_approxim[[3]][,3],color=ID),alpha=alpha_stream, linetype=2)+
      geom_line(aes(x=DtoU,y=four_approxim[[3]][,3],color=ID, alpha=ID), 
                size=1.3,stat="smooth",method = "lm")+
      scale_alpha_manual(values = c(rep(0.3,6),0.9))+
      scale_fill_CUNILLERA(palette = "LGTBI")+
      scale_color_CUNILLERA(palette = "LGTBI")+
      labs(subtitle = "B1 - STcon")+
      xlab("HOBO relative position Upstream-Downstream")+
      ylab("ST Connecticity")+
      scale_y_continuous(limits = c(0,11.5))+
      theme_classic()+theme(legend.position="none"),
    top=""),
  
  #B1 NMDS plot
  arrangeGrob(
    ggplot(Un_NonW_ST_MatrixOut)+
      geom_point(aes(x=x, y=y, fill=ID,size=DtoU), shape=21,colour="grey10", alpha=alpha_stream)+
      stat_ellipse(aes(x=x, y=y, colour=ID),type = "t", linetype=2,size=1, alpha=0.5)+
      scale_fill_CUNILLERA(palette = "LGTBI")+
      scale_color_CUNILLERA(palette = "LGTBI")+
      labs(subtitle = "B1 - STconmat")+
      scale_y_continuous(label=NULL)+
      scale_x_continuous(label=NULL)+
      theme_classic()+theme(legend.position="none",axis.ticks = element_blank())+
      facet_grid(. ~ ID),top = ""),
  
  #B2 line plot  
  arrangeGrob(
    ggplot(data = four_approxim[[4]][,])+
      geom_point(aes(x=DtoU,y=four_approxim[[4]][,3]),color="red",alpha=alpha_HOBOS3, shape=16, size=5,)+
      geom_point(aes(x=DtoU,y=four_approxim[[4]][,3],fill=ID),color="grey30", alpha=alpha_stream, shape=21, size=2)+
      geom_line(aes(x=DtoU,y=four_approxim[[4]][,3],color=ID),alpha=alpha_stream, linetype=2)+
      geom_line(aes(x=DtoU,y=four_approxim[[4]][,3],color=ID, alpha=ID), 
                size=1.3,stat="smooth",method = "lm")+
      scale_alpha_manual(values = c(rep(0.3,6),0.9))+
      scale_fill_CUNILLERA(palette = "LGTBI")+
      scale_color_CUNILLERA(palette = "LGTBI")+
      labs(subtitle = "B2 - STcon")+
      xlab("HOBO relative position Upstream-Downstream")+
      ylab("ST Connecticity")+
      scale_y_continuous(limits = c(0,52000))+
      theme_classic()+theme(legend.position="none"),
    top=""),
  
  #B1 NMDS plot
  arrangeGrob(
    ggplot(Un_WEIG_ST_MatrixOut)+
      geom_point(aes(x=x, y=y, fill=ID,size=DtoU), shape=21,colour="grey10", alpha=alpha_stream)+
      stat_ellipse(aes(x=x, y=y, colour=ID),type = "t", linetype=2,size=1, alpha=0.5)+
      scale_fill_CUNILLERA(palette = "LGTBI")+
      scale_color_CUNILLERA(palette = "LGTBI")+
      scale_y_continuous(label=NULL)+
      scale_x_continuous(label=NULL)+
      labs(subtitle = "B2 - STconmat")+
      theme_classic()+theme(legend.position="none",axis.ticks = element_blank())+
      facet_grid(. ~ ID),top = ""),
  

    get_legend(ggplot(NonW_ST_directed_out)+
                 geom_point(aes(x=DtoU,y=four_approxim[[1]][,3],fill=ID),shape=21, size=5)+
                 scale_fill_CUNILLERA(palette = "LGTBI", name="Stream ID")+
                 theme_classic()+theme(legend.direction = "horizontal",legend.box="vertical")),
    nrow=5 ,ncol=2,widths=c(0.6,1.1))

Charging and depurating biological data

### IMPORTANT: TO PROPERLY RUN THIS SCRIPT YOU NEED TO RUN the following scripts: 

## SpaTemp_HOBOS_treatment.R - The most important one as is the one calculating the ST indices and 
# selecting and creating the HOBOS dataset of interest and then preparing all the data in order to be analysed. 

## Old_HOBOS_calculation.R - Not that important but interesting. Here the Old HOBOS indices are calculated. This 
# indices are calcualted based on the amounts of "1" or "0" of the HOBOS table. They are interesting as they offer
# an extra perspective on the whole ST dynamics. 
## IMPORTANT: If you do not want to use them you will need to remove the parts where they are "called". Generally,
# they are named "Old_HOBOS" or similar names (TotDur, TotNum and TotLeng).

## SpaTemp_comparison.R - It is not essential to run this script. However, it is interesting to check how the 
# HOBOS obtained indices behave among them. Specially to understand what they mean and what they do not mean. So, 
# my recomendation is to run it. 

# SUGGESTED ORDER TO RUN the scripts: 
### 1 SpaTemp_HOBOS_treatment.R
### 2 Old_HOBOS_calculation.R
### 3 SpaTemp_comparison.R
### 4 SpaTemp_Biolcomparison.R

# Cunillera_palette - Available at https://github.com/Cunillera-Montcusi/Cunillera_Pallete
source("C:/Users/David CM/Dropbox/DAVID DOC/LLAM al DIA/CUNILLERA_palette.R")
setwd("C:/Users/David CM/Dropbox/DAVID DOC/LLAM al DIA/1. FEHM coses al DIA/4. Mecodispers Spa-Tem/MECODISPER SpaTem")

# Autumn dataset - November 2019
BDD <- read.csv2("BiolData/Matriz_otoño.csv", sep=";")
# We correct an error as there are two SA4 for October sampling, we sum them. 
BDD[which(BDD$Code=="SA4")[1],3:ncol(BDD)] <- apply(BDD[which(BDD$Code=="SA4"),3:ncol(BDD)],2,sum)
BDD <- BDD[-which(BDD$Code=="SA4")[2],]
BDD <- mutate(BDD, 
              Riera = case_when(
                str_detect(Code,"CA") ~ "C", 
                str_detect(Code,"H") ~ "VH",
                str_detect(Code,"MU") ~ "M",
                str_detect(Code,"R") ~ "R",
                str_detect(Code,"SA") ~ "SA",
                str_detect(Code,"SC") ~ "SC",
                str_detect(Code,"T") ~ "T",
                TRUE ~ "ERROR" )) 

Traits_val <- read.csv2("BiolData/traits.csv", row.names = 1)
bloc<-read.table("BiolData/trait.blocks.txt",h=T)

# TRAITS implementation based on rel. abundances of each trait ####
# Filter species that are present or not present existing 
out <- c()
for (i in 1:length(rownames(Traits_val))) {
  temp_out <- which(rownames(Traits_val)==colnames(BDD%>%dplyr::select(-Date,-Riera,-Code,))[i])
  out[i] <- ifelse(length(temp_out)>0,temp_out,0)
}
Trait_matr_prep<-prep.fuzzy.var(Traits_val[out,],bloc$bloc)
apply(Trait_matr_prep,1,sum)

BDD_croped <- BDD%>%dplyr::select(matches(rownames(Trait_matr_prep)))
rownames(BDD_croped) <-BDD$Code

for(i in 1:nrow(BDD_croped)){BDD_croped[i,]<- BDD_croped[i,]/sum(BDD_croped[i,])}
apply(BDD_croped,1,sum)

Trait_abundances<-as.matrix(BDD_croped)%*%as.matrix(Trait_matr_prep)
  
Trait_abund_f4 <- data.frame("Code"=BDD$Code,"Date"=BDD$Date,"Riera"=BDD$Riera, "f4"=as.data.frame(Trait_abundances)$f4)
Trait_abund_f1 <- data.frame("Code"=BDD$Code,"Date"=BDD$Date,"Riera"=BDD$Riera, "f1"=as.data.frame(Trait_abundances)$f1)

# TRAITS by filtering ####
# Traits addition by  filtering only species corresponding to the traits
Traits_val_edited <- data.frame("taxon"=row.names(Traits_val),Traits_val)

# Active dispersers -- f4==3 and 2
Trait_disp <- Traits_val_edited%>%dplyr::select(taxon,f4)%>%filter(f4==3)%>%
  rows_insert(Traits_val_edited%>%dplyr::select(taxon,f4)%>%filter(f4==2))
BDD_activ <- BDD%>%dplyr::select(Code, Date, Riera,matches(Trait_disp$taxon))

# Passive dispersers -- f1==3 and 2
Trait_disp <- Traits_val_edited%>%dplyr::select(taxon,f1)%>%filter(f1==3)%>%
  rows_insert(Traits_val_edited%>%dplyr::select(taxon,f1)%>%filter(f1==2))
BDD_pass <- BDD%>%dplyr::select(Code, Date, Riera,matches(Trait_disp$taxon))

Active dispersers community metrics calculated

Data_season <- BDD_activ %>%
                dplyr::select(-Date,-Riera)%>%
                pivot_longer(-Code)

colnames(Data_season) <- c("Code","Species","Abund")

#INDIVIDUAL diversity values active ________________________________________________________________________####
by_group_data_season <- group_by(Data_season, 
                                  Code, Species)

data_G_gamma <- by_group_data_season%>%
                group_by(Species)%>%
                summarise(n=sum(Abund))%>%
                mutate(s=1)%>%
                summarise(rich=sum(s))

data_A_alpha <- by_group_data_season%>%
                group_by(Code)%>%
                mutate(s=ifelse(Abund>0,1,0))%>%
                summarise(rich=sum(s))

library(vegan)
table_shannon <- by_group_data_season%>%
                 spread(key =Species,value = Abund, fill=0)%>%
                 bind_cols("Riera"=BDD_activ$Riera)

shan <- diversity(table_shannon%>%ungroup()%>%dplyr::select(-Code,-Riera),index = "shannon")

#SIMMILARITIES active ________________________________________________________________________####

table_shannon_PA <- by_group_data_season%>%
  mutate(PA=ifelse(Abund>0,1,0))%>%
  dplyr::select(-Abund)%>%
  spread(key =Species,value = PA, fill=0)%>%
  bind_cols("Riera"=BDD_activ$Riera)

Bray.distance <- list()
Jac.distance <- list()
for (riera in 1:length(unique(table_shannon$Riera))) {
  names_riera <- unique(table_shannon$Riera[order(table_shannon$Riera)])
  
  # Bray Curtis active
  Bray.distance[[riera]] <- vegdist(
                          decostand(table_shannon%>%
                                      filter(Riera==names_riera[riera])%>%ungroup()%>%dplyr::select(-Code,-Riera),
                          method = "hellinger"),
                          method = "bray")
  
  # Jaccard active
  Jac.distance[[riera]]<-vegdist(
                        decostand(table_shannon_PA%>%
                                    filter(Riera==names_riera[riera])%>%ungroup()%>%dplyr::select(-Code,-Riera),
                        method = "pa"),
                        method = "jaccard")
  }

names(Bray.distance) <- unique(table_shannon$Riera[order(table_shannon$Riera)])
names(Jac.distance) <- unique(table_shannon$Riera[order(table_shannon$Riera)])

BID_output_IND_act <- data.frame(Code=unique(by_group_data_season$Code),rich_f4=data_A_alpha$rich,shan_f4=shan,f4_abun=Trait_abund_f4$f4)
BID_output_DIST_act <- list(Bray.distance,Jac.distance)

Passive dispersers community metrics calculated

Data_season <- BDD_pass %>%
  dplyr::select(-Date,-Riera)%>%
  pivot_longer(-Code)

colnames(Data_season) <- c("Code","Species","Abund")

#INDIVIDUAL diversity values passive ________________________________________________________________________####
by_group_data_season <- group_by(Data_season, 
                                 Code, Species)

data_G_gamma <- by_group_data_season%>%
  group_by(Species)%>%
  summarise(n=sum(Abund))%>%
  mutate(s=1)%>%
  summarise(rich=sum(s))

data_A_alpha <- by_group_data_season%>%
  group_by(Code)%>%
  mutate(s=ifelse(Abund>0,1,0))%>%
  summarise(rich=sum(s))

library(vegan)
table_shannon <- by_group_data_season%>%
  spread(key =Species,value = Abund, fill=0)%>%
  bind_cols("Riera"=BDD_pass$Riera)

shan <- diversity(table_shannon%>%ungroup()%>%dplyr::select(-Code,-Riera),index = "shannon")

#SIMMILARITIES passive ________________________________________________________________________####
table_shannon_PA <- by_group_data_season%>%
  mutate(PA=ifelse(Abund>0,1,0))%>%
  dplyr::select(-Abund)%>%
  spread(key =Species,value = PA, fill=0)%>%
  bind_cols("Riera"=BDD_pass$Riera)

Bray.distance <- list()
Jac.distance <- list()
for (riera in 1:length(unique(table_shannon$Riera))) {
  names_riera <- unique(table_shannon$Riera[order(table_shannon$Riera)])
  
  # Bray Curtis
  Bray.distance[[riera]] <- vegdist(
    decostand(table_shannon%>%
                filter(Riera==names_riera[riera])%>%ungroup()%>%dplyr::select(-Code,-Riera),
              method = "hellinger"),
    method = "bray")
  
  # Jaccard
  Jac.distance[[riera]]<-vegdist(
    decostand(table_shannon_PA%>%
                filter(Riera==names_riera[riera])%>%ungroup()%>%dplyr::select(-Code,-Riera),
              method = "pa"),
    method = "jaccard")
}

names(Bray.distance) <- unique(table_shannon$Riera[order(table_shannon$Riera)])
names(Jac.distance) <- unique(table_shannon$Riera[order(table_shannon$Riera)])

BID_output_IND_pas <- data.frame(Code=unique(by_group_data_season$Code),rich_f1=data_A_alpha$rich,shan_f1=shan,f1_abun=Trait_abund_f1$f1)
BID_output_DIST_pas <- list(Bray.distance,Jac.distance)

Charging geographic and matching biological samples and data logger localization

#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#
# CHARGING GEOGRAPHIC COORDINATES & MATCHING SAMPLES AND HOBOS ####
#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#

# We charge the field Samples and filter them according to the ones actually taken (from BDD)  ####
SampSites <- read.table("BiolData/Sites_1.txt",sep = ",",header = T)
out <- c()
for (match in 1:length(BDD$Code)) {
 a <- which(SampSites$Site==BDD$Code[match])
 out[match] <- ifelse(length(a)==0,0,a)
}
SampSites <- SampSites[out,]

# We call the ST indices to MATCH them  ####
## You need to run the previous script in order to built the HOBOS dataset (a list with the HOBOS info for each river)
Sites_list # They are charged from SpaTemp_HOBOS_treatment.R
HOBOS_sites# They are charged from SpaTemp_HOBOS_treatment.R

# This two other objects are also built in the "SpaTemp_HOBOS_treatment.R" script. They contain the names 
## of each river and the order upstream to downstream within it. 
# They are used to Identify and locate river position between HOBOS and real Samples
HOB_riv_ID<- Stream_order$ï..Riera
ups_dos <- Stream_order$UtoD

# Calculation of OLD HOBOS values to include them in the comparisons. 
source("Old_HOBOS_calculation.R")
Old_HOBOS_comb

# Geographical matching ####
# Converting Sampling SITES to apropriate geog. data 
str(SampSites)
a <- c()
b <- c()
for (posit in 1:nrow(SampSites)) {
a[posit] <-as.numeric(sub("(.{2})(.*)", "\\1.\\2", SampSites$lat[posit]))
b[posit] <-as.numeric(sub("(.{1})(.*)", "\\1.\\2", SampSites$long[posit]))
}
SampSites$lat<- a
SampSites$long<- b

# Converting HOBOS SITES to appropriate geog. data 
### *IMPORTANT STEP ####
#Here is where we include the calculated individual ST indices by "left_join"
Sites_list_comb <- rbind(Sites_list[[1]],Sites_list[[2]],Sites_list[[3]],Sites_list[[4]],Sites_list[[5]],
                         Sites_list[[6]],Sites_list[[7]] )
Sites_list_comb <- cbind(Sites_list_comb, "ID"=HOB_riv_ID,"DtoU"=ups_dos)%>%
                   mutate(ID_UpDo=paste(ID,"_",DtoU,sep = ""))
Sites_list_comb <- Sites_list_comb%>%
                   left_join(NonW_ST_directed_out)%>%
                   left_join(WEIG_ST_directed_out)%>%
                   left_join(Un_NonW_ST_out)%>%
                   left_join(Un_WEIG_ST_out)%>%
                   left_join(Old_HOBOS_comb)
a <- c()
b <- c()
for (posit in 1:nrow(Sites_list_comb)) {
  a[posit] <-as.numeric(sub("(.{2})(.*)", "\\1.\\2", as.character(round(Sites_list_comb$Latitud[posit]))))
  b[posit] <-as.numeric(sub("(.{1})(.*)", "\\1.\\2", as.character(round(Sites_list_comb$Longitud[posit]))))
}
Sites_list_comb$Latitud<- a
Sites_list_comb$Longitud<- b

# Plot to check that the values match. 
par(mfrow=c(1,2))
plot(Sites_list_comb$Longitud, Sites_list_comb$Latitud)
plot(SampSites$long,SampSites$lat)
par(mfrow=c(1,1))

# Match the values
library(RANN) 
Matching_HOBOS <- nn2(data=Sites_list_comb[,3:4], query = SampSites[,2:3], k=1, radius = 1)[[1]] 

### *MANUAL CHECKING ####
# Here you must make sure that the assignation is correct between Samp_ID, ID, Riera, and DtoU
HOB_BDD_match <- as.data.frame(cbind(Sites_list_comb[Matching_HOBOS,], "Samp_ID"=SampSites$Site)) # Check assignation
data.frame(
HOB_BDD_match$Riera,
HOB_BDD_match$Codi_HOBO,
HOB_BDD_match$DtoU,
HOB_BDD_match$Samp_ID)

# Order the sites by upstream downstream directionality
HOB_BDD_match <- HOB_BDD_match%>%arrange(Riera,DtoU)

#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#
# MERGING ALL THE VALUES TOGETHER ####
#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#

#Individual
HOB_BDD_match <- HOB_BDD_match%>%left_join(BID_output_IND, by=c("Samp_ID"="Code"))

#Matrices
### *IMPORTANT STEP ####
#Here is where we include the calculated matrix ST indices
NonW_Dir <- list()
WEIG_Dir<- list()
Un_NonW<- list()
Un_WEIG<- list()
for (riera in 1:length(unique(NonW_ST_matrix_out_out))) {
  names_riera <- unique(Sites_list_comb$ID)
  
  values_places <- unlist(HOB_BDD_match%>%filter(ID==names_riera[riera])%>%dplyr::select(DtoU))
  
  NonW_Dir[[riera]] <- as.dist(t(NonW_ST_matrix_out_out[[riera]][values_places,values_places]))
  WEIG_Dir[[riera]] <- as.dist(t(WEIG_ST_matrix_out_out[[riera]][values_places,values_places]))
  Un_NonW[[riera]] <- as.dist(t(Un_NonW_ST_matrix_out_out[[riera]][values_places,values_places]))
  Un_WEIG[[riera]] <- as.dist(t(Un_WEIG_ST_matrix_out_out[[riera]][values_places,values_places]))
}

### *MANUAL CHECKING ####
# We eliminate R and G matrices due to the lack of sites for this two sampling places
## For diversities (Bray and Jaccard) eliminate positions 2 and 4
## For ST matrices eliminate positions 3
HOB_BDD_match_matrix <- list("BrayCurtis_Act"=BID_output_DIST$Active[[1]][-c(2,4)],
                             "Jaccard_Act"=BID_output_DIST$Active[[2]][-c(2,4)],
                             "BrayCurtis_Pas"=BID_output_DIST$Passive[[1]][-c(2,4)],
                             "Jaccard_Pas"=BID_output_DIST$Passive[[2]][-c(2,4)],
                             "NonW_Dir"=NonW_Dir[-c(3)], 
                             "WEIG_Dir"=WEIG_Dir[-c(3)], 
                             "Un_NonW"=Un_NonW[-c(3)], 
                             "Un_WEIG"=Un_WEIG[-c(3)])

# We eliminate R and G in Individual due to the lack of samples from this period
HOB_BDD_match <- HOB_BDD_match%>%filter(Samp_ID !=c("R4","G1"))

# FINAL MERGED VALUES - INDIVIDUAL AND DISTANCE (MATRIX) BASED 
HOB_BDD_match
HOB_BDD_match_matrix

### *SMALL EDITIONS TO ELIMINATE Ocl/Acl/Bet ####
HOB_BDD_match <- HOB_BDD_match%>%dplyr::select(-c(NonW_Dir_Ocl,NonW_Dir_Acl,NonW_Dir_Bet,
                                           WEIG_Dir_Ocl,WEIG_Dir_Acl,WEIG_Dir_Bet,
                                           Un_NonW_Ocl,Un_NonW_Acl,Un_NonW_Bet,
                                           Un_WEIG_Ocl,Un_WEIG_Acl,Un_WEIG_Bet))

Richness

plots_HOB_BDD_total <- list()
plots_HOB_BDD <- list()
model_HOB_BDD_results <- list()
for (col_var in 1:ncol(HOB_BDD_match%>%dplyr::select(-Riera, -Codi_HOBO,
                                              -Latitud,-Longitud,
                                              -ID, -DtoU,-ID_UpDo,
                                              -Samp_ID, -TotDur, -TotNum, -TotLeng, 
                                              -rich_f4, -shan_f4,-f4_abun,
                                              -rich_f1,-shan_f1,-f1_abun))){

  
  variable_x_temp <- HOB_BDD_match%>%dplyr::select(-Riera, -Codi_HOBO,
                                       -Latitud,-Longitud,
                                       -ID, -DtoU,-ID_UpDo,
                                       -Samp_ID, -TotDur, -TotNum, -TotLeng, 
                                       -rich_f4, -shan_f4,-f4_abun,
                                       -rich_f1,-shan_f1,-f1_abun)

  variable_x_name <- colnames(variable_x_temp)[col_var]
  variable_x <- variable_x_temp[,col_var]
  
  variable_y_temp <- HOB_BDD_match%>%dplyr::select(rich_f4,rich_f1)
  plot_counter <- c(col_var, col_var+ncol(variable_x_temp))
  
for (variable in 1:ncol(variable_y_temp)) {
    
  variable_y <- variable_y_temp[,variable]
  variable_y_model <- unlist(variable_y)
  variable_y_name <- colnames(variable_y_temp)
  
  factors <- HOB_BDD_match%>%dplyr::select(ID, DtoU,ID_UpDo)
  Id_random <- unlist(factors$ID)
  
  model <- lme(log(variable_y_model+1)~log(variable_x+1), random = ~1|Id_random)
  results <- round(as.numeric(c(summary(model)[[4]]$fixed, summary(model)[[20]][2,5])),2)
  #Results table
  model_table <- round(summary(model)[[20]],3)
  rownames(model_table) <- c("Intercept", "STconmat")
  model_table<- rownames_to_column(as.data.frame(model_table))
  model_table<- cbind(Scenario=names_scenarios[col_var],Dispersion=Dispersal_group[variable],model_table)
  model_HOB_BDD_results[[plot_counter[variable]]] <- model_table 
  
  dataset <- cbind(factors, "X_var"=variable_x, variable_y)
  dataset$pred_values <- predict(model)   
  
plots_HOB_BDD[[plot_counter[variable]]] <- ggplot(dataset,aes(x=log(X_var+1), y=log(variable_y+1),fill=ID,colour=ID))+
    geom_point(color="grey30", alpha=0.5, shape=21, size=2)+
    geom_smooth(method = "lm", colour="grey60",alpha=0.1, se = TRUE, linetype=2)+
    geom_abline(slope =results[2],intercept = results[1], colour="black", size=2)+
    scale_fill_manual(values =  CunilleraPAL_corrected)+
    scale_color_manual(values =  CunilleraPAL_corrected)+
    xlab(paste("log(STcon)"))+
    ylab(paste("log(Richness)"))+
    labs(caption = paste("Intercept=",results[1],"Slope=",results[2],"p-value=",results[3]))+
    labs(subtitle=paste(names_scenarios[col_var],"vs",variable_y_name[variable]))+
    theme_classic()+
    theme(legend.position="none")
  }
}

model_result <- rbind(model_HOB_BDD_results[[1]],model_HOB_BDD_results[[2]],
                      model_HOB_BDD_results[[3]],model_HOB_BDD_results[[4]],
                      model_HOB_BDD_results[[5]],model_HOB_BDD_results[[6]],
                      model_HOB_BDD_results[[7]],model_HOB_BDD_results[[8]])
colnames(model_result)[[3]] <- "Fixed effects"
write.table(model_result, "Table_Results/Biol_Rich_STcon.txt", sep = ",")

legend_plots<- get_legend(ggplot(dataset)+
                            geom_point(aes(x=X_var,y=variable_y ,fill=ID),shape=21, size=6)+
                            scale_fill_manual(values = CunilleraPAL_corrected, name="Stream ID")+
                            theme_classic()+theme(legend.direction = "vertical",legend.box="vertical"))

grid.arrange(
  plots_HOB_BDD[[1]],plots_HOB_BDD[[2]],plots_HOB_BDD[[3]],plots_HOB_BDD[[4]],
  legend_plots,
  plots_HOB_BDD[[5]],plots_HOB_BDD[[6]], plots_HOB_BDD[[7]],plots_HOB_BDD[[8]],
  nrow=2,ncol=5, top="Richness")

plots_HOB_BDD_total[[1]] <- plots_HOB_BDD

Shannon-Wiener

plots_HOB_BDD <- list()
model_HOB_BDD_results <- list()
for (col_var in 1:ncol(HOB_BDD_match%>%dplyr::select(-Riera, -Codi_HOBO,
                                              -Latitud,-Longitud,
                                              -ID, -DtoU,-ID_UpDo,
                                              -Samp_ID, -TotDur, -TotNum, -TotLeng, 
                                              -rich_f4, -shan_f4,-f4_abun,
                                              -rich_f1,-shan_f1,-f1_abun))){
  
  
  variable_x_temp <- HOB_BDD_match%>%dplyr::select(-Riera, -Codi_HOBO,
                                            -Latitud,-Longitud,
                                            -ID, -DtoU,-ID_UpDo,
                                            -Samp_ID, -TotDur, -TotNum, -TotLeng, 
                                            -rich_f4, -shan_f4,-f4_abun,
                                            -rich_f1,-shan_f1,-f1_abun)
  
  variable_x_name <- colnames(variable_x_temp)[col_var]
  variable_x <- variable_x_temp[,col_var]
  
  variable_y_temp <- HOB_BDD_match%>%dplyr::select(shan_f4,shan_f1)
  plot_counter <- c(col_var, col_var+ncol(variable_x_temp))
  
  for (variable in 1:ncol(variable_y_temp)) {
    
    variable_y <- variable_y_temp[,variable]
    variable_y_model <- unlist(variable_y)
    variable_y_name <- colnames(variable_y_temp)
    
    factors <- HOB_BDD_match%>%dplyr::select(ID, DtoU,ID_UpDo)
    Id_random <- unlist(factors$ID)
    
    model <- lme(log(variable_y_model+1)~log(variable_x+1), random = ~1|Id_random)
    results <- round(as.numeric(c(summary(model)[[4]]$fixed, summary(model)[[20]][2,5])),2)
    #Results table
    model_table <- round(summary(model)[[20]],3)
    rownames(model_table) <- c("Intercept", "STconmat")
    model_table<- rownames_to_column(as.data.frame(model_table))
    model_table<- cbind(Scenario=names_scenarios[col_var],Dispersion=Dispersal_group[variable],model_table)
    model_HOB_BDD_results[[plot_counter[variable]]] <- model_table 
    
    dataset <- cbind(factors, "X_var"=variable_x, variable_y)
    dataset$pred_values <- predict(model)   
    
    plots_HOB_BDD[[plot_counter[variable]]] <- ggplot(dataset,aes(x=log(X_var+1), y=log(variable_y+1),fill=ID,colour=ID))+
      geom_point(color="grey30", alpha=0.5, shape=21, size=2)+
      geom_smooth(method = "lm", colour="grey60",alpha=0.1, se = TRUE, linetype=2)+
      geom_abline(slope =results[2],intercept = results[1], colour="black", size=2)+
      scale_fill_manual(values =  CunilleraPAL_corrected)+
      scale_color_manual(values =  CunilleraPAL_corrected)+
      xlab(paste("log(STcon)"))+
      ylab(paste("log(Shannon)"))+
      labs(caption = paste("Intercept=",results[1],"Slope=",results[2],"p-value=",results[3]))+
      labs(subtitle=paste(names_scenarios[col_var],"vs",variable_y_name[variable]))+
      theme_classic()+
      theme(legend.position="none")
  }
}

model_result <- rbind(model_HOB_BDD_results[[1]],model_HOB_BDD_results[[2]],
                      model_HOB_BDD_results[[3]],model_HOB_BDD_results[[4]],
                      model_HOB_BDD_results[[5]],model_HOB_BDD_results[[6]],
                      model_HOB_BDD_results[[7]],model_HOB_BDD_results[[8]])
colnames(model_result)[[3]] <- "Fixed effects"
write.table(model_result, "Table_Results/Biol_ShaWei_STcon.txt", sep = ",")

legend_plots<- get_legend(ggplot(dataset)+
                            geom_point(aes(x=X_var,y=variable_y ,fill=ID),shape=21, size=6)+
                            scale_fill_manual(values = CunilleraPAL_corrected, name="Stream ID")+
                            theme_classic()+theme(legend.direction = "vertical",legend.box="vertical"))

grid.arrange(
  plots_HOB_BDD[[1]],plots_HOB_BDD[[2]],plots_HOB_BDD[[3]],plots_HOB_BDD[[4]],
  legend_plots,
  plots_HOB_BDD[[5]],plots_HOB_BDD[[6]], plots_HOB_BDD[[7]],plots_HOB_BDD[[8]],
  nrow=2,ncol=5, top="Shannon")

plots_HOB_BDD_total[[2]] <- plots_HOB_BDD

Trait abundance

plots_HOB_BDD <- list()
model_HOB_BDD_results <- list()
for (col_var in 1:ncol(HOB_BDD_match%>%dplyr::select(-Riera, -Codi_HOBO,
                                              -Latitud,-Longitud,
                                              -ID, -DtoU,-ID_UpDo,
                                              -Samp_ID, -TotDur, -TotNum, -TotLeng, 
                                              -rich_f4, -shan_f4,-f4_abun,
                                              -rich_f1,-shan_f1,-f1_abun))){
  
  
  variable_x_temp <- HOB_BDD_match%>%dplyr::select(-Riera, -Codi_HOBO,
                                            -Latitud,-Longitud,
                                            -ID, -DtoU,-ID_UpDo,
                                            -Samp_ID, -TotDur, -TotNum, -TotLeng, 
                                            -rich_f4, -shan_f4,-f4_abun,
                                            -rich_f1,-shan_f1,-f1_abun)
  
  variable_x_name <- colnames(variable_x_temp)[col_var]
  variable_x <- variable_x_temp[,col_var]
  
  variable_y_temp <- HOB_BDD_match%>%dplyr::select(f4_abun,f1_abun)
  plot_counter <- c(col_var, col_var+ncol(variable_x_temp))
  
  for (variable in 1:ncol(variable_y_temp)) {
    
    variable_y <- variable_y_temp[,variable]
    variable_y_model <- unlist(variable_y)
    variable_y_name <- colnames(variable_y_temp)
    
    factors <- HOB_BDD_match%>%dplyr::select(ID, DtoU,ID_UpDo)
    Id_random <- unlist(factors$ID)
    
    model <- lme(log(variable_y_model+1)~log(variable_x+1), random = ~1|Id_random)
    results <- round(as.numeric(c(summary(model)[[4]]$fixed, summary(model)[[20]][2,5])),2)
    #Results table
    model_table <- round(summary(model)[[20]],3)
    rownames(model_table) <- c("Intercept", "STconmat")
    model_table<- rownames_to_column(as.data.frame(model_table))
    model_table<- cbind(Scenario=names_scenarios[col_var],Dispersion=Dispersal_group[variable],model_table)
    model_HOB_BDD_results[[plot_counter[variable]]] <- model_table 
    
    dataset <- cbind(factors, "X_var"=variable_x, variable_y)
    dataset$pred_values <- predict(model)   
    
    plots_HOB_BDD[[plot_counter[variable]]] <- ggplot(dataset,aes(x=log(X_var+1), y=log(variable_y+1),fill=ID,colour=ID))+
      geom_point(color="grey30", alpha=0.5, shape=21, size=2)+
      geom_smooth(method = "lm", colour="grey60",alpha=0.1, se = TRUE, linetype=2)+
      geom_abline(slope =results[2],intercept = results[1], colour="black", size=2)+
      scale_fill_manual(values =  CunilleraPAL_corrected)+
      scale_color_manual(values =  CunilleraPAL_corrected)+
      xlab(paste("log(STcon)"))+
      ylab(paste("log(Trait abundance)"))+
      labs(caption = paste("Intercept=",results[1],"Slope=",results[2],"p-value=",results[3]))+
      labs(subtitle=paste(names_scenarios[col_var],"vs",variable_y_name[variable]))+
      theme_classic()+
      theme(legend.position="none")
  }
}

model_result <- rbind(model_HOB_BDD_results[[1]],model_HOB_BDD_results[[2]],
                      model_HOB_BDD_results[[3]],model_HOB_BDD_results[[4]],
                      model_HOB_BDD_results[[5]],model_HOB_BDD_results[[6]],
                      model_HOB_BDD_results[[7]],model_HOB_BDD_results[[8]])
colnames(model_result)[[3]] <- "Fixed effects"
write.table(model_result, "Table_Results/Biol_TraitAb_STcon.txt", sep = ",")

legend_plots<- get_legend(ggplot(dataset)+
                            geom_point(aes(x=X_var,y=variable_y ,fill=ID),shape=21, size=6)+
                            scale_fill_manual(values = CunilleraPAL_corrected, name="Stream ID")+
                            theme_classic()+theme(legend.direction = "vertical",legend.box="vertical"))

grid.arrange(
  plots_HOB_BDD[[1]],plots_HOB_BDD[[2]],plots_HOB_BDD[[3]],plots_HOB_BDD[[4]],
  legend_plots,
  plots_HOB_BDD[[5]],plots_HOB_BDD[[6]], plots_HOB_BDD[[7]],plots_HOB_BDD[[8]],
  nrow=2,ncol=5, top="Trait abundance")

plots_HOB_BDD_total[[3]] <- plots_HOB_BDD

Matching Bray & Jaccard with STconmat values

output <- matrix(nrow=length(unlist(HOB_BDD_match_matrix$BrayCurtis_Act)))
for (lis_elem in 1:length(HOB_BDD_match_matrix)) {
  out <- c()
  for (river in 1:length(HOB_BDD_match_matrix$BrayCurtis_Act)) {
    col_names_river <- names(HOB_BDD_match_matrix$BrayCurtis_Act)
    
    out_temp <- c()
    out_temp <- cbind(as.matrix(HOB_BDD_match_matrix[[lis_elem]][[river]])[
      upper.tri(as.matrix(HOB_BDD_match_matrix[[lis_elem]][[river]]))])
    out_temp <- as.data.frame(cbind("ID"=rep(col_names_river[river],nrow(out_temp)),out_temp))
    out<- rbind(out, out_temp)
  }
  colnames(out)[2] <- names(HOB_BDD_match_matrix)[lis_elem]
  output[,1] <- out[,1]
  output <- as.data.frame(cbind(output, out))
  output <- output%>%dplyr::select(-ID)
  output[,lis_elem+1] <- as.numeric(output[,lis_elem+1])
}
colnames(output)[1] <- "ID"

STmatrix_BiolDissim <- output

Bray-Curtis simmilarity

plots_HOB_BDD <- list()
model_HOB_BDD_results <- list()
for (col_var in 1:ncol(STmatrix_BiolDissim%>%dplyr::select(-ID,
                                                    -BrayCurtis_Act, -Jaccard_Act,
                                                    -BrayCurtis_Pas, -Jaccard_Pas))){
  
  
  variable_x_temp <- STmatrix_BiolDissim%>%dplyr::select(-ID,
                                             -BrayCurtis_Act, -Jaccard_Act,
                                             -BrayCurtis_Pas, -Jaccard_Pas)

  variable_x_name <- colnames(variable_x_temp)[col_var]
  variable_x <- variable_x_temp[,col_var]
  
  variable_y_temp <- STmatrix_BiolDissim%>%dplyr::select(BrayCurtis_Act,BrayCurtis_Pas)
  plot_counter <- c(col_var, col_var+ncol(variable_x_temp))
  
  for (variable in 1:ncol(variable_y_temp)) {
    
    variable_y <- variable_y_temp[,variable]
    variable_y_model <- unlist(variable_y)
    variable_y_name <- colnames(variable_y_temp)
    
    factors <- STmatrix_BiolDissim%>%dplyr::select(ID)
    Id_random <- unlist(factors$ID)
    
    model <- lme(variable_y_model~log(variable_x+1), random = ~1|Id_random)
    results <- round(as.numeric(c(summary(model)[[4]]$fixed, summary(model)[[20]][2,5])),2)
    #Results table
    model_table <- round(summary(model)[[20]],3)
    rownames(model_table) <- c("Intercept", "STconmat")
    model_table<- rownames_to_column(as.data.frame(model_table))
    model_table<- cbind(Scenario=names_scenarios[col_var],Dispersion=Dispersal_group[variable],model_table)
    model_HOB_BDD_results[[plot_counter[variable]]] <- model_table 
    dataset <- cbind(factors, "X_var"=variable_x, variable_y)
    dataset$pred_values <- predict(model)   
    
    plots_HOB_BDD[[plot_counter[variable]]] <- ggplot(dataset,aes(x=log(X_var+1), y=variable_y,fill=ID,colour=ID))+
      geom_point(color="grey30", alpha=0.5, shape=21, size=2)+
      geom_smooth(method = "lm", colour="grey60",alpha=0.1, se = TRUE, linetype=2)+
      geom_abline(slope =results[2],intercept = results[1], colour="black", size=2)+
      scale_fill_manual(values =  CunilleraPAL_corrected)+
      scale_color_manual(values =  CunilleraPAL_corrected)+
      xlab(paste("log(STconmat)"))+
      ylab(paste("Bray curtis"))+
      labs(caption = paste("Intercept=",results[1],"Slope=",results[2],"p-value=",results[3]))+
      labs(subtitle=paste(names_scenarios[col_var],"vs",variable_y_name[variable]))+
      theme_classic()+
      theme(legend.position="none")
  }
}

model_result <- rbind(model_HOB_BDD_results[[1]],model_HOB_BDD_results[[2]],
                      model_HOB_BDD_results[[3]],model_HOB_BDD_results[[4]],
                      model_HOB_BDD_results[[5]],model_HOB_BDD_results[[6]],
                      model_HOB_BDD_results[[7]],model_HOB_BDD_results[[8]])
colnames(model_result)[[3]] <- "Fixed effects"
write.table(model_result, "Table_Results/Biol_BrayC_STconmat.txt", sep = ",")

legend_plots<- get_legend(ggplot(dataset)+
                            geom_point(aes(x=X_var,y=variable_y,fill=ID),shape=21, size=6)+
                            scale_fill_manual(values = CunilleraPAL_corrected, name="Stream ID")+
                            theme_classic()+theme(legend.direction = "vertical",legend.box="vertical"))

grid.arrange(
  plots_HOB_BDD[[1]],plots_HOB_BDD[[2]],plots_HOB_BDD[[3]],plots_HOB_BDD[[4]],
  legend_plots,
  plots_HOB_BDD[[5]],plots_HOB_BDD[[6]], plots_HOB_BDD[[7]],plots_HOB_BDD[[8]],
  nrow=2,ncol=5, top="Bray Crutis")

plots_HOB_BDD_total[[4]] <- plots_HOB_BDD

Jaccard

plots_HOB_BDD <- list()
model_HOB_BDD_results <- list()
for (col_var in 1:ncol(STmatrix_BiolDissim%>%dplyr::select(-ID,
                                                    -BrayCurtis_Act, -Jaccard_Act,
                                                    -BrayCurtis_Pas, -Jaccard_Pas))){
  
  variable_x_temp <- STmatrix_BiolDissim%>%dplyr::select(-ID,
                                                  -BrayCurtis_Act, -Jaccard_Act,
                                                  -BrayCurtis_Pas, -Jaccard_Pas)
  
  variable_x_name <- colnames(variable_x_temp)[col_var]
  variable_x <- variable_x_temp[,col_var]
  
  variable_y_temp <- STmatrix_BiolDissim%>%dplyr::select(Jaccard_Act,Jaccard_Pas)
  plot_counter <- c(col_var, col_var+ncol(variable_x_temp))
  
  Dispersal_group <- c("Active", "Passive")
  for (variable in 1:ncol(variable_y_temp)) {
    
    variable_y <- variable_y_temp[,variable]
    variable_y_model <- unlist(variable_y)
    variable_y_name <- colnames(variable_y_temp)
    
    factors <- STmatrix_BiolDissim%>%dplyr::select(ID)
    Id_random <- unlist(factors$ID)
    
    model <- lme(variable_y_model~log(variable_x+1), random = ~1|Id_random)
    results <- round(as.numeric(c(summary(model)[[4]]$fixed, summary(model)[[20]][2,5])),2)
    #Results table
    model_table <- round(summary(model)[[20]],3)
    rownames(model_table) <- c("Intercept", "STconmat")
    model_table<- rownames_to_column(as.data.frame(model_table))
    model_table<- cbind(Scenario=names_scenarios[col_var],Dispersion=Dispersal_group[variable],model_table)
    model_HOB_BDD_results[[plot_counter[variable]]] <- model_table 
    
    dataset <- cbind(factors, "X_var"=variable_x, variable_y)
    dataset$pred_values <- predict(model)   
    
    plots_HOB_BDD[[plot_counter[variable]]] <- ggplot(dataset,aes(x=log(X_var+1), y=variable_y,fill=ID,colour=ID))+
      geom_point(color="grey30", alpha=0.5, shape=21, size=2)+
      geom_smooth(method = "lm", colour="grey60",alpha=0.1, se = TRUE, linetype=2)+
      geom_abline(slope =results[2],intercept = results[1], colour="black", size=2)+
      #geom_smooth(method = "lm",alpha=0.2, se = T, fill="grey50", colour="black", size=2)+
      scale_fill_manual(values =  CunilleraPAL_corrected)+
      scale_color_manual(values =  CunilleraPAL_corrected)+
      xlab(paste("log(STconmat)"))+
      ylab(paste("Jaccard"))+
      labs(caption = paste("Intercept=",results[1],"Slope=",results[2],"p-value=",results[3]))+
      labs(subtitle=paste(names_scenarios[col_var],"vs",variable_y_name[variable]))+
      theme_classic()+
      theme(legend.position="none")
    #facet_grid(.~ID)  
  }
}

model_result <- rbind(model_HOB_BDD_results[[1]],model_HOB_BDD_results[[2]],
                      model_HOB_BDD_results[[3]],model_HOB_BDD_results[[4]],
                      model_HOB_BDD_results[[5]],model_HOB_BDD_results[[6]],
                      model_HOB_BDD_results[[7]],model_HOB_BDD_results[[8]])
colnames(model_result)[[3]] <- "Fixed effects"
write.table(model_result, "Table_Results/Biol_Jacc_STconmat.txt", sep = ",")


legend_plots<- get_legend(ggplot(dataset)+
                            geom_point(aes(x=X_var,y=variable_y,fill=ID),shape=21, size=6)+
                            scale_fill_manual(values = CunilleraPAL_corrected, name="Stream ID")+
                            theme_classic()+theme(legend.direction = "vertical",legend.box="vertical"))

grid.arrange(
  plots_HOB_BDD[[1]],plots_HOB_BDD[[2]],plots_HOB_BDD[[3]],plots_HOB_BDD[[4]],
  legend_plots,
  plots_HOB_BDD[[5]],plots_HOB_BDD[[6]], plots_HOB_BDD[[7]],plots_HOB_BDD[[8]],
  nrow=2,ncol=5, top="Jaccard")

plots_HOB_BDD_total[[5]] <- plots_HOB_BDD